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ABSTRACT 

Burg's  algorithm  for  Maximum  Entropy  autoregressive  spectral  estimation  is 
analyzed  for  the  cases  of  one  and  two  complex  sinusoidal  signals  in  additive  white 
noise.  For  the  latter  case  are  found  two  biases  which  can  account  for  the  line 
splitting  and  line  shifting  that  occur  in  simulation  studies  when  the  SNR  is  very 
high.  These  biases  vanish  completely  if  the  two  complex  sinusoids  are  in  phase 
quadrature  at  the  middle  of  the  data  record;  If  there  is  an  integral  number  of 
half-cycles  of  difference  frequency  contained  in  the  data  record,  then  the  power 
level  of  the  spectral  estimate  will  be  biased  although  the  effects  believed  to  cause 
splitting  and  shifting  will  be  eliminated.  Results  of  simulation  studies  to  support 
these  conjectures  are  presented. 

f~ - — 


1.  INTRODUCTION 

The  Burg  algorithm  for  the  autoregressive  spectral  analysis  of  time- 
series  data1*2*3,  sometimes  referred  to  as  the  maximum-entropy  method  (MEM), 
is  known  to  be  inappropriate  for  the  case  of  sinusoidal  signals  in  additive 
white  noise.  This  inappropriateness  had  been  demonstrated  both  theoretically 1,-6 
and  in  practice7-9.  A  theoretically  correct  model1*-8  for  the  generating  process 
for  an  N-pole  complex  sinusoidal  signal  in  additive  white  noise  is  an  N-pole, 
N-zero  network  with  Identical  gain  weights  in  its  feedback  (pole)  and  feed¬ 
forward  (zero)  parts,  being  excited  with  a  white-noise  input  (see  Figure  1). 
Autogresslve  analysis  models  the  generating  process  for  the  data  as  an  all-pole 
network  excited  by  white  noise.  Since  a  zero  in  the  generating  process 
network  can  be  simulated  exactly  only  by  an  infinite  number  of  poles,  it  is 
clear  that  when  autoregressive  analysis  is  used,  in  principle  an  infinite  set 
of  either  autocorrelation  or  time  series  data  must  be  used  in  order  to  achieve 
correct  results. 


r(2))  (y(i) 


Figure  1.  Network  for  generating  N  complex  sine  waves  in  additive  white  noise  from  complex  white  noise  input 

(after  Figure  4. 1  of  Reference  61 


Fougere’  has  stated  that,  in  the  high  signal-to-noise  ratio  (SNR)  case, 
the  Burg  algorithm  is  overconstrained.  In  the  case  of  simulation  studies,  this 
overconstraint  causes  errors  in  the  estimated  frequencies  of  spectral  lines 
and  the  false  splitting  of  spectral  lines  known  to  have  been  generated  by  a 
single  pole  (or  a  pair  of  poles,  in  the  case  of  real  signals).  Fougere  has 
developed  an  algorithm  which  avoids  this  phenomenon,  but  the  algorithm  is  based 
on  a  gradient-search  technique  which  lacks  the  intrinsic  efficiency  of  the 
unmodified  Burg  algorithm. 

In  spite  of  its  known  limitations,  the  Burg  algorithm  is  often  used 
because  of  its  computational  efficiency.  This  report  explores  analytically 
the  expected  response  of  the  Burg  algorithm  to  time-series  data  comprising  one 
or  two  complex  sinusoids,  with  and  without  the  presence  of  additive  white 
noise.  It  is  shown  that  only  in  very  special  cases  does  the  Burg  algorithm 
lead  to  the  same  results  as  are  achieved  when  the  true  autocorrelation  functions 
of  the  signals  are  known. 


REVIEW  OF  AUTOREGRESSIVE  SPECTRAL  ANALYSIS 


Autoregressive  spectral  analysis  is  based  on  the  idea  that.  If  it  is 
somehow  possible  to  design  a  feedforward  (all-zero)  filter  which  has  as  its 
input  the  data  to  be  analyzed,  and  has  as  its  output  random  white  noise,  then 
the  power  spectrum  of  the  input  data  is  given  by  the  reciprocal  of  the  power 
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transfer  function  of  the  filter.  Since  this  filter  accounts  for  all  the 
predictability  Inherent  in  the  input  signal  and  has  as  its  output  only  unpre¬ 
dictable  random  white  noise,  it  is  often  referred  to  as  a  prediction-error 
filter  (PEF) . 

There  are  several  well-known  techniques  for  estimating  the  PEF 
corresponding  to  a  given  set  of  data.  When  only  amplitude  time-series  data 
are  available,  most  of  these  depend  upon  estimates  of  the  autocorrelation 
function  derived  from  the  time-series  data.  The  Burg  algorithm,  however, 
attempts  to  avoid  possible  biases  or  inconsistencies  in  such  estimates  of  the 
autocorrelation  function  by  deriving  an  estimate  of  the  PEF  coefficients 
directly  from  the  data. 

in  this  report  the  algorithm  for  generating  the  I’EF  when  the  true 
autocorrelation  function  is  known  is  reviewed  in  Section  2.1  and  the  Burg 
algorithm  is  reviewed  in  Section  2.2.  The  results  from  these  sections  are 
then  used  to  generate  the  sets  of  PEFs  corresponding  to  the  cases  of  one  and 
two  complex  sinusoids  in  additive  white  noise,  and  the  properties  of  these 
sets  of  PEFs  are  compared  and  contrasted. 

2.1  THE  KNOWN-AUTOCORRELATION  (KA)  ALGORITHM1 0,11 

Let  it  be  assumed  that  N  equlspaced  samples  R(n)  of  the  complex  auto¬ 
correlation  function  have  been  given,  for  n*=0,l,2 . N-l,  where  N  may  be 

finite  or  infinite.  It  is  assumed  that  the  Nyquist  sampling  criterion  has 
been  met.  Then  the  system  of  equations  to  be  solved  is 


M 

-  E  R(k-m)a(m,M)  =■  P(M)6(k)  (1) 

m*=0 


0  <_  k  <  M 
0  <  M  <  N-l 


where  the  a(m,M)s  for  m=>0,l,2 . M  are  sets  of  PEF  coefficients,  each  value 

for  M  denoting  a  different  set;  the  P(M)s  are  called  the  output  error  powers 
and  are  real;  and  6(k)  =  1  if  k-0,  6(k)  =  0  otherwise,  is  the  Kroenecker  delta 
function.  In  order  to  maintain  proper  scaling,  the  leading  terms  of  the  PEFs, 
a(0,M),  M*0,1, . . . ,N-1,  are  set  equal  to  -1  by  definition. 

Since  negative  indices  for  R(k-m)  occur  in  the  set  of  equations  (1), 
it  is  necessary  to  note  that  R(-n)  *  R*(n),  where  the  asterisk  (*)  denotes 
complex  conjugation.  This  last  fact  allows  (1)  to  be  rewritten  in  the  alter¬ 
native  form: 

M 

-  E  R(k,+m)ot*(m,M)  -  P(M)6(k')  (2) 

m*0 

-M  1  k*  1  0 
0  <  M  <  N-l 
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Equations  (2)  imply  that  the  same  result  is  obtained  if  the  complex  conjugate 
of  a  PEF  is  applied  to  the  time-reversed  autocorrelation  data.  This  "reverse- 
conjugate"  symmetry  is  used  in  the  derivation  of  the  Burg  algorithm. 

If  the  set  of  linear  equations  (1)  for  a  particular  value  of  M  is 
written  in  matrix  form,  it  can  be  seen  that  the  matrix  of  autocorrelation 
samples  [R(k,m)],  where  R(k,m)  =  R(k-m),  is  an  MxM  Toeplitz  matrix.  Therefore 
(1)  can  be  solved  by  applying  the  algorithm  developed  by  Levinson.  Robinson 
and  Durbin,  otherwise  known  as  the  Levinson  recursion.  Following11  the 
recursive  solution  to  (1)  can  be  written  as: 


P(0)  =  R(0) 

(3a) 

M-l 

a(M,M) 

«  -  E  a (m ,M-1 ) R (M-m) /P (M-l ) 

(3b) 

m=0 

a(m,M) 

=  a(m,M-l)  -  a(M,M)a*(M-m  ,M-1) 

(3c) 

P(M) 

-(1  -  |a(M,M)|2)P(M-l) 

(3d) 

m  *  1,2,... ,M-1 
M  «=  1,2,... 


Note  that  at  each  successive  stage  of  the  recursion,  the  introduction  of  one 
new  autocorrelation  sample,  R(M),  generates  but  one  Independent  value  a(M,M) 
for  the  M^-order  PEF;  all  other  coefficients  of  the  M^-order  PEF  are 
determined  from  linear  combinations  of  the  coefficients  of  the  (M-l)t^l-order 
PEF  and  their  complex  conjugates,  using  a(M,M)  as  indicated  by  (3c). 

The  ct(M,M)s  are  sometimes  referred  to  as  the  reflection  coefficients, 
because  of  the  analogy  of  their  appearance  in  (3d)  with  a  similar  equation 
which  occurs  in  the  theory  of  a  signal  propagating  through  a  layered  medium 
and  being  partially  reflected  at  each  layer  interface. 

In  executing  the  recursion  of  equations  (3a-d)  it  can  occur  (at  least 
in  theory)  that,  for  some  particular  value  of  M,  say  MQ,  P(M0)  *  0.  This 
implies  that  |a(M0,M0) |  ■  1.  This  condition  can  arise  only  in  the  case  where 
the  signal  being  analyzed  can  be  modelled  as  M0  complex  sinusoids  with  no 
additive  noise  (see  Sections  3.1  and  3.2);  in  general  M0  is  not  finite.  In 
particular,  M0  cannot  be  finite  when  additive  white  noise  is  present4. 

For  each  order  M  of  PEF,  an  estimate  of  the  power  spectrum  based  on 
(M+l)  values  of  the  autocorrelation  function  is  given  by 


XKa«o,m) 


_ _ PM _ 

|M 

E  a(m,M)exp(-jmu) 
m-0 


(A) 
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where  oj  is  the  normalized  angular  frequency  in  radians  with  -ir  <  u  <_  ir,  and 
the  subscript  KA  refers  to  the  known-autocorrelation  case.  As  discussed  in 
Section  2.0  and  as  examination  of  (1)  will  indicate,  the  PEFs  are  "spiking" 
or  whitening  filters,  since  all  but  one  of  their  output  values  are  zero.  The 
power  spectrum  of  such  an  output  signal  is  independent  of  frequency,  or  "white". 
The  denominator  in  the  right-hand  term  of  (4)  is  the  power  transfer  function 
of  the  PEF,  which  if  multiplied  by  Xk^(id,M),  the  estimated  spectrum  of  the 
signal,  yields  the  constant,  white-noise  spectrum  P(M).  Thus,  since  both  P(M) 
and  the  power  spectrum  of  the  PEF  can  be  calculated,  the  signal  power  spectrum 
Xka(u),M)  can  be  estimated  from  (4)  for  successively  higher  orders  of  PEF. 


2.2  THE  BURG  ALGORITHM 

The  Burg  algorithm  is  a  procedure  for  estimating  the  reflection 
coefficients  directly  from  a  set  of  time-series  amplitude  data.  It  avoids 
the  biases  introduced  into  the  spectral  estimate  when  the  autocorrelation 
is  estimated  from  the  data  and  the  known-autocorrelation  algorithm  is  then 
applied;  however,  as  is  shown  in  Sections  4.1  and  4.2,  the  Burg  algorithm 
introduces  biases  of  its  own  sort. 

Let  it  be  assumed  that  a  set  of  N  time-series  amplitude  data  x(n), 
n=0,l ,2, . . . ,N-1  have  been  given,  and  that  <x(n)>  =  0,  where  the  brackets  <  > 
denote  expected  value.  The  PEFs  are  derived  sequentially.  Each  successively 
higher  order  PEF  is  applied  to  the  data  in  both  directions  simultaneously,  and 
the  average  of  its  forward  and  backward  output  error  powers  is  minimized  by 
adjusting  only  its  reflection  coefficient.  The  remaining  coefficients  of 
each  PEF  depend  on  the  sequence  of  reflection  coefficients  through  the 
functional  relationship  defined  by  (3c).  The  motivation  for  this  procedure 
is  its  analogy  with  that  defined  by  (1),  (2),  and  (3a-d)  . 

Following  e.g.,3,11-1<*  and  taking  proper  note  of  the  occurrences  of 
complex  conjugation  in  the  complex  data  case,  the  Burg  algorithm  can  be 
written  in  a  lattice-filter  formulation: 


and 


fM(n) 


M-l 


Vn) 


M-l 


f  (n) 
o 


(n)  - 

B(M,M)bM_i(n-1) 

(5a) 

(n-1) 

-  6*(M,M)fM_1(n) 

(5b) 

n  =  M ,M+1 , . 

. . ,N-1 

M  -  1,2,... 

,N-1 

*  x(n) 

-bo(n) 

(5c) 

n  =  0,1,2, . 

. ,,N-1 

(See  Figure  2,  where  z  ^  denotes  the  unit  time-delay  operator.)  The  series 

fu(n)  is  the  output  from  the  M^-order  PEF  applied  to  the  input  signal  x(n) 
M 


6 


in  Che  forward  direction,  and  is  expressed  in  terms  of  the  output  series  from 
the  (M-l)t^-stage  of  the  lattice.  The  series  b^(n)  is  the  output  from  the 
Mth_order  PEF,  conjugated  and  applied  to  the  input  data  in  the  reverse  or 
backward  direction,  and  again  is  expressed  in  terms  of  the  output  series  from 
the  (M-l)t*1-stage  of  the  lattice.  The  B(M,M)s  are  the  reflection  coefficients. 


fo(n)  f|(n)  fg(n)  fM(n) 


b0(n)  b,(n)  b2(n)  bM(n) 


Figure  2.  Basic  all-zero  lattice  network 
(after  Reference  11) 


The  sum  of  the  forward  and  backward  output  error  energies  at  each 
stage  of  the  lattice  is  given  by 


N-l 

E(M)  -  1  (|fM(n)|2+  I bM<n)  | 2)  (6) 

n=m 

M  =  0,1,2 . N-l 

The  formula  for  computing  B(M,M),  the  Burg  estimate  of  the  Mtbl-order  reflection 
coefficient,  is  derived  by  substituting  eqns.  (5a)  and  (5b)  into  (6),  setting 
[ 3E (M) / 3 B* (M ,M) ]  =  0  and  solving  for  B(M,M). 


2  X  bVl(#-1)fM-lW 

B(M,M)  ■=  — - 2^ -  (7) 

Z  <|bM_lCn-1)|2  +  |fM-l<n>|2> 


Notice  the  similarity  of  (7)  to  a  single-lag  unwlndowed  cross-correlation  of 
the  forward  and  backward  output  series. 


7 


In  order  to  obtain  spectral  estimates.  It  Is  usual  to  let  the  output 
error  powers  11  be  defined  as 


w(0)  -  E(0)/(2N)  (8a) 

and 


tt(M)  =  (1  -  |  $(M,M)  |  2)tt(M-1)  (8b) 

M  =  1,2, . . . ,N-1 

by  analogy  with  (3a)  and  (3d)  respectively.  Then,  letting  the  P(M,M)s  be 
defined  by 


fl(m,M)  =  fl(m.M-l)  -  P(M ,M) p* (M-m ,M-1)  (9) 

m  ■  1,2, .. . ,M-1 

by  analogy  with  (3c),  and  P(0,M)  =  -1  by  definition,  the  Mth-order  Burg  power 
spectrum  estimate  XgOd.M)  is  given  by 


Xb(io,N) 


tt(M) _ 

li  p" 

£  $(m,M)exp(-jmuj) 
m=0 


by  analogy  with  (4) . 


(10) 


2.3  THE  ZEROES  OF  THE  PREDICTION-ERROR  FILTERS  (PEFs) 

By  employing  standard  z-transform  techniques,  the  z-transforms  of  the 
PEFs  for  the  KA  case  and  the  Burg  case  can  be  written  as  polynomials  in  the 
complex  variable  z.  For  the  KA  case,  this  polynomial  is  F  (z,M)  where 

KA 

M 

Fka(z,M)  *=  1  -  £  a(m,M)z"in  (11) 

m*l 

Then  (4)  can  be  rewritten  as 

1^(0), M)  -  P(M)/tFKA(exp(ja))>M]|2  (12) 


where  the  denominator  is  the  squared  magnitude  of  FKA(z,M)  evaluated  around 
the  unit  circle  (|z|  =  1).  Similarly,  for  the  Burg  case,  the  polynomial  is 
Fb(z,M),  where 


r 


H 


M 

F_(z,M)  -  1  -  Z  3(m,M)z"in 


(13) 


and  (10)  can  be  rewritten  as 


Xb(gu,M)  =  Tr(M)/|FB[exp(jaj),M]|2 


(14) 


It  Is  apparent  from  (12)  and  (14)  that  if  any  zeroes  of  FKA(z,M)  or 
I’b(z.M)  lie  near  or  on  the  unit  circle,  then  the  magnitude  of  the  spectral 
estimators  X|^(w,M)  or  Xp(w,M)  will  be  large  at  locations  on  the  unit  circle 
in  the  vicinity  of  such  zeroes.  Conversely,  zeroes  lying  close  to  the  origin 
of  the  complex  z-plane  will  have  little  influence  on  the  peaks  of  the  spectrum, 
but  will  affect  its  magnitude  away  from  the  peaks.  Thus  some  insight  into  the 
character  of  an  autoregressive  spectral  estimate  can  be  gained  by  studying  the 
locations  of  the  zeroes  of  its  associated  PEF. 


3.  THE  ONE-POLE  COMPLEX  SINUSOID  CASE 

The  formula  for  a  complex  signal  (xjfn))  consisting  of  a  single 
complex  sinusoid  in  the  presence  of  additive  white  noise  is  given  by 


x^n)  =  exp(jnw1)  +  e(n)  (15) 

where  n  is  any  positive  or  negative  integer,  or  zero;  A1  is  the  complex 
amplitude  of  the  complex  sinusoid;  is  the  angular  frequency  of  the  complex 
sinusoid,  normalized  so  that  -tt  <  <_  it;  and  e(n)  is  additive  white  noise 

having  the  property 


<e*(n)e(n+k)>  =  |e|26(k) 


(16) 


where  |e|2  is  used  to  denote  the  variance  of  e(n)  and  6(k)  is  again  the 
Kroenecker  delta  function. 


3.1  THE  KA  ESTIMATE  OF  THE  PEF 

The  autocorrelation  function  R^k)  of  the  signal  defined  by  (15)  is 

given  by 


Rx(k)  *  |AjJ 2  exp(jku>1)  +  | e| 2<S(k) 


(17) 
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where  k  is  any  positive  or  negative  integer,  or  zero,  and  in  general  R(k)  is 
defined  by 


R(k)  *=  <x*(n)x(n+k)>  (18) 

Substitution  of  (17)  into  (3a)  gives  the  result 

Pj (0)  =  | e| 2 (oj  +  1)  (19) 

where  oj  =  |  A-^  |  2  / 1  e  |  2  is  the  signal-to-noise  ratio  (SNR).  If  there  is  no 
noise,  i.e.  |ej2  =  0,  then  Pi(0)  =  | | 2 • 

The  following  general  result  can  be  derived  from  (3b),  (3c)  and  (3d) 
when  | e | 2  t  0  and  |A^|?  +  0: 


oij  (m,M)  =  [a2/(Ma2+l)]exp(jmu)1) 

(20a) 

=  1 1/  (M+o-^)  ]  exp  ( jmtOj) 

(20b) 

P1(M)  =  | e | 2  i[(M+l)aJ+l]/[MoJ+l]| 

(21) 

M  =  1,2,3, 
m  =  1,2,3 . M 

In  the  noise-free  case  (|e|2  -*■  0  and  aj  -*•“>) ,  (20a)  or  (20b)  imply 
that  a]_(l,l)  =  exptjwj)  so  that  lajd.l)!2  =  1.  Then  (3d)  implies  that 
Pj^(l)  =  0.  Thus,  for  this  special  case,  the  sequences  defined  by  (20)  and  (21) 
terminate  at  M=l.  Otherwise  for  this  signal  model  the  sequences  are  infinite, 
with  Pi (M)  -*■  |e|2(l+M-l)  as  Mcrf  ~. 

The  z-transform  of  the  KA  PEF  (cf.  (11)) 


F^z.M)  =  1  -  [ai/(MoJ+l)]  Z  z"*  expOmiOj)  (22) 

m»l 


where  the  superscript  1  denotes  a  one-pole  complex  sinusoidal  signal.  For 
M=l,  F^)(z,l)  has  a  zero  at  zQ“[af/(a2+l)  ]exp(ju1) .  For  M-2,  F^)(z,2)  has 
zeroes  at . 


z 

o 


a2  ±  a 


/9a2  +  4 
1  1 


exp  (JUj) 


(23) 
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so  that  for  high  SNR  (1  <<  o|  <  “) ,  the  zeroes  occur  at  approximately  1.0 
exp(jwj)  and  -0.5  exp(ju)j),  and  for  low  SNR  (o|  <<  1),  the  zeroes  occur  at 
approximately  +°i  expfjto^). 

Some  algebraic  manipulation  shows  that,  for  the  product  MoJ  suffi¬ 
ciently  greater  than  1,  an  approximate  root  zQ  of  (22)  is  given  by 


z 

o 


(M-l)  Ma?  I  exp 


f] 


(j^) 


(24) 


This  means  that,  as  MoJ  -*■  °°,  the  estimated  location  of  the  pole  corresponding 
to  the  single  complex  sinusoid  asymptotically  approaches  the  correct  location 
exp(jaj^)  on  the  complex  z-plane  along  a  radius  oriented  at  an  angle  correspond¬ 
ing  to  the  true  frequency  of  the  signal.  This  is  true  even  when  of  «  1. 

Numerical  solutions  of  (24)  show  that  the  other  zeroes  tend  to 
distribute  themselves  with  approximately  uniform  angular  separation  and 
approximately  constant  radius  inside  the  unit  circle  so  as  to  account  for  the 
uniform  spectrum  of  the  additive  white  noise.  The  radius  at  which  the  zeroes 
occur  varies  inversely  with  the  SNR. 


3.2  THE  BURG  ESTIMATE  OF  THE  PEF 

Let  it  be  assumed  that  N  samples  of  the  signal  defined  by  (15)  have 
been  given: 


x1(n)  =  A1  exp(jnoi^)  +  e(n) 

n  =  0,1,2, . . . ,N-1 


(25) 


Substituting  (25)  into  (6)  by  using  (5c)  and  then  taking  the  expected  value 
of  E^(0)  yields 


<E1(0)>  =  2N | e | 2(a£+l)  (26) 

It  is  not  so  easy  to  calculate  <3^(1, 1)>,  since  examination  of  (7) 
shows  that  it  is  necessary  to  derive  the  expected  value  of  a  quotient  of 
correlated  random  variables.  In  general,  this  requires  that  the  statistics 
of  the  random  variables  be  specified  and  the  problem  be  solved  numerically. 
This  will  not  be  done  here. 

However,  for  sufficiently  high  SNR  (e.g.,  >  100  or  20  dB)  the 

approximation  (l+x)“l  ~  1-x,  where  x  ~  e(n)/A^,  can  be  used  t_  approximate 
the  denominator  of  (7).  Then  the  following  approximaf  isult,  which  is 
independent  of  the  statistical  distribution  of  the  white  noise  e(n),  is 
obtained : 


|l/[l+B1(l,N)o1~2]|  exp(jw1) 


(27) 


<B1(1,1)>  * 


where 


B1(1,N)  -  (Nz-2) / (N-l) 2  (28) 

and  2  B^d.N)  >  1  for  2  <_  N  <  °°.  Comparison  of  (27)  with  (20b)  for  M=1 
shows  that  for  high  SNR  Bj_(l,l)  Is  a  biased  estimate  of  ai(l,l).  However, 
Pl(l,l)  correctly  estimates  the  angular  frequency  wj  of  the  single  complex 
sinusoid,  and  the  bias  term  B^(1,N)  monotonically  approaches  unity  as  N 
becomes  large. 

From  (26) ,  (27)  and  (8a)  and  (8b)  it  can  be  shown  that 

<tt1(0)>  =  |  e  |  2  (o£+l)  (29) 

and,  to  the  same  degree  of  approximation  as  was  used  to  obtain  (27) 

<it1(1>>  -  2B1(l,N)|e|2  (30) 

Comparison  of  (30)  with  (21)  shows  that  tt^(1)  is  a  biased  estimator  of  Pj(l), 
since  for  oj  large,  Pj.(l)  ~  2 1 e J  2 . 

The  result  (30)  Implies  that  the  lattice-filter  outputs  fj^(n)  and 
b^(n)  as  defined  by  (5a)  and  (5b)  have  low  SNR,  since  their  expected  power 
is  at  most  only  a  factor  of  four  greater  than  the  additive  white  noise  power 
| G | 2 .  Therefore  it  would  again  be  necessary  to  specify  the  detailed  statistics 
of  the  additive  noise  before  the  higher  order  reflection  coefficients  could  be 
estimated.  This  will  not  be  done  here. 


3.3  DISCUSSION  OF  THE  ONE-POLE  COMPLEX  SINUSOID  CASE 

The  results  of  Sections  3.1  and  3.2  show  that,  for  high  SNR  (o|  >  100) 
the  first  order  (M*l)  PEF  generated  by  the  Burg  algorithm  is  biased  as  compared 
to  the  PEF  generated  by  the  KA  technique.  This  bias,  however,  monotonically 
decreases  as  the  number  of  data  N  is  increased,  and  both  the  KA  and  the  Burg 
algorithms  correctly  estimate  the  frequency  of  the  complex  sinusoid. 

It  is  impossible  to  investigate  the  properties  of  the  Burg  PEFs  for 
the  low  SNR  case,  or  for  orders  higher  than  M»1  for  the  high  SNR  case,  without 
specifying  the  statistical  distribution  of  the  additive  white  noise.  This 
problem  has  not  been  considered  here. 


4.  THE  TWO-POLE  COMPLEX  SINUSOID  CASE 

The  sampled  signal  X2(n)  consisting  of  two  complex  sinusoids  in  the 
presence  of  additive  white  noise  is  given  by 
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x2(n)  *  Ax  exp(jnw^)  +  A2  exp(jnu>2)  +  e(n)  (31) 

til 

where  =  |  A^|  exp(j4>ic)  Is  the  complex  amplitude  of  the  k  sinusoid,  ^  is 
its  arbitrary  initial  phase  at  n=0,  and  cu^  is  its  angular  frequency,  normalized 
so  that  -it  <  0)^  <_  tt,  for  k=l,2.  e(n)  is  additive  white  noise,  as  in  (15).  It 
is  apparent  that  if  A2=Aj  and  0)2  =  -0)j_  then  X2(n)  is  a  sampled  real  sinusoid 
in  additive  white  (complex)  noise. 

Equation  (31)  can  be  written  in  a  form  which  will  be  subsequently 
more  tractable: 


x9(n)  =  A  exp [  j  (nu>  +$  ) ] 
c  o  o  o 

X  jr  exp[  j  (nAurt-A<}>)]  +  r  1  exp[-j  (nAaH-A<J>)  ]|  +  e(n)  (32) 

where  A0  **  I A^2 | ^  is  the  geometric  mean  of  the  magnitudes  of  the  two 
amplitudes;  r  =  [ |AjJ / | A2 | l*5  is  the  square  root  of  the  ratio  of  the  magnitudes 
of  the  amplitudes;  <J>0  =  <4>^-Kj>2) / 2  is  the  arithmetic  mean  of  the  initial  phases; 
A<f>  =  (4)^— 4>2) / 2  is  one-half  the  difference  between  the  two  phases;  o)Q  *  (oj^+o)2)/2 
is  the  arithmetic  mean  of  the  two  angular  frequencies;  and  Aco  ■  (oi^— 0*2) / 2  is 
one-half  the  difference  between  the  two  frequencies. 


4.1  THE  KA  ESTIMATE  OF  THE  PEF 

The  autocorrelation  function  R2(k)  of  the  signal  defined  by  (32)  is 

given  by 

R2(k)  =  A2(r2+r  2)exp(jkuio)  [coskAw  +  jp(r)sinkAw)  +  |  e |  2 6(k)  (33) 

where 

p(r)  =  (r2-r-2)/ (r2+r~2)  (34) 

and  k  is  any  positive  or  negative  integer,  or  zero.  Substitution  of  (33)  into 
(3a)  gives  the  result 

P2(0)  =  |e|2(o2+1)  (35) 

_2 

where  o|  *  Ao(r2+r  )/|e|2  is  the  SNR.  If  there  is  no  noise,  then  P2(0)  = 
Ao(r2+r~2)  is  the  signal  power. 

From  (33),  (35),  (3b)  and  (3d),  the  following  general  results  for 

the  reflection  coefficient  and  the  output  error  power  for  M®1  can  be  derived 

when  I  £  I  2  ^  0  and  A  ^  0 : 

1  '  o 
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a2(l,l)  -  exp(jwQ) j[cosAw  +  jp(r)sinAw] /[I  +  a~2]}  (36) 

and 

P2(l)  =  A2(r2+r“2){[l-p2(r)]sin2Aw  +  o~2(2-k>22)}/{l-Kr22}  (37) 

For  high  SNR  (o|  »  1),  the  first  reflection  coefficient  is  given 
approximately  by 


a2(l,l)  =  exp(jWQ)  [cosAw  +  jp(r)sinAw]  (38) 

For  r=l,  p(r)  =  0  and  the  zero  of  F^)(z,l)  (see  Section  2.3)  is  located  at 
zG  =  cosAw  exp(jwD),  which  lies  on  a  radius  oriented  at  the  mean  angular 
frequency  wG.  This  zero  moves  towards  the  limiting  values  of  exp(jwj)  as  r-**» 
and  p(r)  1,  or  expQti^)  as  and  p(r)  -+•  -1,  and  the  single  complex  sinusoid 
case  is  approached  in  each  case. 

The  output  error  power  is,  using  (38),  (35)  and  (8b), 

P2(l)  =  (2Aq  sinAw)2/ (r2+r-2)  (39) 

which  is  essentially  the  signal  power  attenuated  by  the  factor  l2sinAw/ (r2+r~  )]2. 
This  factor  is  unity  when  r=l  and  Aw-+tt/2,  and  decreases  as  r  or  Aw  deviate 
from  these  values. 

For  low  SNR  (o|  <<  1),  the  reflection  coefficient  and  the  output 
error  power  are  given  by 

a2(l,l)  *  exp(jwo)o2  (cosAw  +  jp(r)sinAw]  (AO) 


and 


P2(l)  -  |e| 2  +  A2 (r2+r~2)  (41) 

Thus  the  frequency  estimated  from  the  location  of  the  zero  of  the  M=1  order 
PEF  lies  in  the  range  bounded  by  wQ  ±  Aw,  and  the  output  error  power  is 
essentially  the  unattenuated  signal  plus  noise  power. 

For  M=2,  (33),  (36),  and  (3b)  can  be  combined  to  yield 


a2(2,2)  =  -exp(j2wo) 


A 

[l-p2(r)]8in2Aw  -  [cos2Aw  +  jp(r)sin2Aw) 


-2, 


[l-p2(r)]sin2Aw  +  a/[2+a,‘] 
and  (37),  (42)  and  (3d)  can  be  combined  to  yield 


(42) 
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2 [  1— p2  (r) ] sin2Aoj(2  +  cos2Aw)  +  o~2(3+o~2) 

P,(2)  =  | e|  2  - 5 - ^ - —  (43) 

f  1-p2  (r)  ]  sin2A<a  +  oy  (2+0^  ) 

For  high  SNR  (cr^  »  1)  and 

o2  [l-p2(r)]sin2Aw  »  |cos2Aw  +  jp(r)sin2Ae»)|  (44) 

(42)  reduces  to 

a2(2,2)  =  -exp(j2too)  (45) 

It  can  be  shown  from  (45),  (36)  and  (3c)  that 

a2(l,2)  =  exp(jaio)  2cosAw  (46) 

so  that,  in  the  limit  as  the  left-hand  side  of  (44)  approaches  infinity, 
the  zeroes  of  Fj^p  approach  z0  =  expf j (wo+Aw) ] ,  the  true  locations  of  the 
poles  of  the  complex  sinusoids.  In  this  same  limit,  | ot2  (2 , 2 1  2  =  1,  so  that 
P2(2)=0  and  the  recursion  terminates. 

For  the  intermediate  case,  where  >>  1  but  (44)  is  not  satisfied, 

i.e., 

o|[l-p2(r))sinAw  «  |cos2Aw  +  jp(r)sin2Aw|  (47) 

(42)  can  be  reduced  to 

a2(2,2)  =  exp(j2oJo)cosAo)[0.5  cosAuj  +  jp(r)sinAu)]  (48) 

and  it  can  be  shown  from  (48) ,  (38)  and  (3c)  that 

a2(l,2)  =  0.5  exp(ju)Q)  [cosAw  +  jp(r)sinAa)]  (49) 

so  that  the  zeroes  of  Fijp(z,2)  occur  at  [cosAw  +  jp(r)sinAo)]exp(j<i)0)  and 
-0.5[cosAw  +  jp(r)sinAo)Jexp(jto0) .  Comparison  of  these  results  with  (23), 
which  gives  the  zero  locations  of  the  PEF  for  a  one-pole  signal  and  M*l, 
shows  that  the  intermediate  case  result  for  the  two-pole  signal  is  very 
similar  to  that  for  the  one-pole  signal  at  high  SNR.  Comparison  with  (38) 
shows  that  one  of  the  zeroes  of  F^)(z,2)  remains  the  same  as  that  of  fX^)(z,1) 
at  high  SNR;  i.e.,  the  two  poles  are  estimated  as  one  by  the  M*2  PEF  if^44) 
is  not  satisfied. 
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For  the  low  SNR  case  (o2  «  1) ,  (42)  reduces  to 

a2(2,2)  =  exp(j2u>o)  jo|[cos2Ati)  +  jp(r)sin2Au)]|  (50) 

and  P2(2)  is  the  same  as  P2(l)  as  given  by  (41)  .  Since  for  this  case  both 
| f»2 (2,2)  |  and  | c»2 ( 1 » 1 )  |  are  proportional  to  o^,  then  from  (3c)  it  is  clear 
that,  to  first  order  in  o^. 


a2(l,2)  =  a2(l,l)  (51) 

where  o^d,  1)  is  given  by  (40).  The  zeroes  of  Ty^  (z,2)  in  this  case  occur 
at  approximately  z  =  02|±[cos2Au)  +  jp(r)sin2Ao)]iS  +  c^tcosAu)  +  jp(r)sinAo)i) /l\ 
exp(j0)o).  For  r=l  and  thus  p(r)=0,  the  zeroes  occur  at  z  =  02  j±[cos2Aio]'5  + 
02[cosAw] /l\ exp(j0)o) ,  which  lie  close  to  the  origin  of  the  complex  z-plane, 
on  a  diameter  of  the  unit  circle  passing  through  the  point  z  =  exp(ju)0).  For 
r  »  1  (p(r)  +  1)  or  r  «  1  (p(r)  -*■  -1)  the  zeroes  tend  to  the  locations 
z  =  C2[±l  +  02/2] exp(ju)i)  or  z  5  021*1  +  02/2] exp(ja>2)  respectively.  Thus  it 
is  apparent  that  for  low  SNR,  the  spectral  estimate  corresponding  to  M=2  is 
incapable  of  resolving  the  spectral  peaks  corresponding  to  the  poles  of  the 
two  complex  sinusoidal  signals. 

There  appears  to  be  no  straightforward  recursion  formula  for  the  KA 
PEF  coefficients,  as  in  the  case  of  the  single  sinusoid  example  of  Section  3.1. 
Therefore,  following  this  approach,  it  is  not  easy  to  determine  the  behavior 
of  the  KA  PEFs  in  the  case  of  large  M  and,  in  particular,  whether  resolution 
of  the  two  sinusoids  is  to  be  expected  for  the  product  M sufficiently  large, 
independent  of  the  value  of  o|.  This  problem,  however,  has  been  solved  using 
powerful  matrix  techniques,  by  Marple6. 


4.2  THE  BURG  ESTIMATE  OF  THE  PEF 

Let  it  again  be  assumed  (cf.  Section  3.2)  that  N  samples  of  the  signal 
defined  by  (32)  have  been  given: 


x2(n)  =  Aq  exp[j(ntoo  +  <t>Q)]{r  exp[  jnAurt-A<|>)  ]  +  r_1exp[-j  (nAuH-A<f>)]}+  e(n)  (52) 

n  =  0,1,2 . N-l 

Substituting  (52)  into  (6)  and  (7),  using  (5c),  and  then  applying  (8a)  and 
taking  expected  values  with  respect  to  the  additive  white  noise  only  yields 


<n2(0)>  •=  A*(r2+r~2)[l  +  2  cosA$mld  G(N,Aw)/(r2+r-2)]  +  |e|2 
for  the  expected  signal-plus-noise  power.  Here 


G(N,Aw) 


sin(NAu)) 
N  sinAui 


(53) 


(54) 
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is  the  common  grating-function  frequency  response  of  a  normalized,  uniformly 
weighted  discrete  Fourier  transform  of  N  data,  and 

A<j)mid  =  (N-l)Aio  +  2A<J>  (55) 

is  the  phase  difference  between  the  two  complex  sinusoidal  components,  reckoned 
at  the  middle  of  the  data  set.  Note  that  there  may  or  may  not  be  a  datum  at 
the  middle  of  the  data  set,  according  to  whether  N  is  an  odd  or  even  integer, 
respectively.  Also  note  that  A<J)  has  not  been  averaged,  but  rather  is  assumed 
to  be  a  fixed  parameter  of  the  particular  set  of  data  being  analyzed.  This 
assumption  corresponds  to  the  usual  practical  case,  where  only  one  set  of  data 
is  available. 

Again,  it  is  not  easy  to  calculate  the  expected  values  of  the 
reflection  coefficients  unless  the  assumption  of  high  SNR  is  made.  In  that 
case  the  same  approximations  can  be  made  as  in  the  derivation  of  (26)  to  get 

cosAw  +  jp(r)sinAw  +  2  cosA<J>  .  G(N-l,Aw)/(r2+r  ) 

<g„(l,l)>  =  exp(ju)  )  - — 75 - 1? - 

1  +  2  cosA<f>  ..cosAo)  G(N-l,Aw)/(r2+r  )  +  cr,  [1  +  B0(1,N)] 

mia  l  l 

(56) 

where  the  bias  term  B2(1,N)  is  of  order  (N-l)  1  and  is  given  by  eqn.  (Al)  of 
Appendix  A. 

Comparison  of  (56)  and  (36)  shows  that  for  high  SNR  Is  a 

biased  estimate  of  (*2(1,1)  unless  N-*50.  For  infinite  SNR  and  finite  N, 
however,  f32(l,l)  becomes  an  unbiased  estimate  of  02(1,1)  if 

cosAcj)  ,  =  0  (57) 

mia 

or 

G (N-l, Aw)  =  0  (58) 

Similarly,  comparison  of  (53)  and  (35)  shows  that  II2(0)  is  a  biased  estimate 
of  P2(0)  unless  either  (57)  is  satisfied  or 

G(N,Aw)  -  0  (59) 

It  is  impossible  to  satisfy  (58)  and  (59)  simultaneously  for  N  finite,  but 
when  (57)  is  satisfied  the  Burg  spectral  estimate  (14)  is  unbiased  for  M*1 
and  infinite  SNR. 

Progressing  now  to  the  second  stage  (M=2)  of  the  Burg  algorithm,  it 
is  found  that  the  algebra  becomes  all  but  intractable  unless  the  condition 
of  infinite  SNR  is  assumed.  For  this  special  case,  (56),  (52),  (5a-c)  and 
(7)  can  be  used  to  derive  the  following  expression  for  32(2,2): 


B2(2,2) 


-exp(j2a)o)  X  ||l  ~  2  cosA4>mjd  G(N-2,Aw)/(r2+r  2)j 

-2  cosA<|>  .  ,G(N-l,Aaj)  X  icosAdi  .  ,G(N-2,Aw)  [cosAui  +  jp(r)sinAu>] 
mla  j  mid  -  \ 

-  2  cosAu)/(r2+r~ 


+  cos2A((>midG2  (N-l,Au))  X  |[cos2Ao)  +  jp(r)sin2Aco] 


-2cosA<|>  .  ,G (N-2 , Aw) /  (r 2+r  2) 
mla 


>! 


-  2cosA<J>m^dG(N-2,Aa))/(r2+r  2)j 

-2cosA<j>mid  cosAwG(N-l,Ad))  X  |cosA<J>midG  (N-2  ,Auj)  -  2/(r2+r  2)| 

+  cos2A<J>  ,  ,G2  (N-l,Au>)  X  '  2cosA<t>  cos2A(jjG(N-2,Aw)/ (r2+r 
mid  I  mid 


Examination  of  (60)  shows  that,  even  for  infinite  SNR,  the  "correct" 
value  of  -exp(j2<i)0)  for  the  reflection  coefficient  <82(2, 2)>  is  not  realized 
for  N  finite  unless  either  (57)  or  (58)  is  satisfied.  Realization  of  either 
of  these  conditions  for  infinite  SNR  will,  as  examination  of  (60)  shows,  cause 
the  magnitude  of  the  reflection  coefficient  to  be  unity.  Thus  it  can  be 
inferred  that,  for  sufficiently  high  SNR,  the  two  crucial  factors  affecting 
the  reflection  coefficient  computed  using  the  Burg  algorithm  are  the  phase 
difference  between  the  two  complex  sinusoids  at  the  middle  of  the  data  record, 
A<f>mid>  an<*  t*ie  number  of  cycles  of  difference  frequency  contained  in  the 
finite  length  data  record. 


4.3  DISCUSSION  OF  THE  TWO-POLE  COMPLEX  SINUSOID  CASE 

The  results  of  Section  4.1  and  4.2  show  that  even  for  infinite  SNR, 
the  first  (M=l)  and  the  second  (M=2)  order  PEFs  generated  by  the  Burg 
algorithm  are  biased  as  compared  to  the  PEFs  generated  by  the  KA  technique. 

It  is  clear  from  the  appearance  of  the  grating  function  (54)  in  the  equations 
(56)  and  (60)  for  the  first  and  second-order  reflection  coefficients  that  the 
magnitude  of  such  biases  will  have  inverse  dependence  on  N,  the  length  of  the 
data  record . 

The  effect  of  this  bias  is  to  reduce  the  magnitude  of  the  reflection 
coefficient  and  thus  to  allow  significant  levels  of  uncancelled  signal  energy 
to  propagate  beyond  the  stage  M*=2  in  the  Burg  algorithm.  Then  PEFs  of 
successively  higher  order  can  be  based  on  this  coherent  "leakage"  signal. 
However,  whenever  one  of  the  criteria  described  by  (57)  or  (58)  is  satisfied, 
no  significant  coherent  leakage  signal  is  propagated  beyond  the  stage  M*2. 

It  is  conjectured  that  it  is  the  presence  or  absence  of  this  coherent  leakage 
signal  beyond  the  stage  M=2  that  determine  whether  or  not  line  splitting  will 
occur  for  PEFs  of  some  higher  order.  Results  both  of  previously  published7*8 
and  new  simulation  studies  support  this  conjecture,  as  indicated  in  Section  5. 
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5.  RESULTS  OF  SOME  SIMULATION  STUDIES 

In  this  Section  are  presented  the  results  of  some  studies  of  the 
performance  of  the  complex  Burg  algorithm  for  the  analysis  of  signals  known 
to  be  comprised  of  two  complex  sinusoids  in  the  presence  of  very  weak  additive 
complex  white  noise  (02  =  77  dB) .  These  studies  parallel  and  extend  a  set  of 
studies  performed  by  Fougere  et  al.8  using  the  real-arithmetic  Burg  algorithm 
to  estimate  spectra  of  a  single  real  sine-wave  signal  in  the  presence  of  very 
weak  additive  real  white  noise. 

Tt  will  he  necessary  to  make  comparisons  between  complex  data  of  the 
form  (52)  and  real  data  xs(n)  of  the  form 

xs(n)  =  sin(n(i)s  +  <j)s>  +  fs(n)  (61) 

n  *  0,1 ... . ,N-1 

which  is  comprised  of  N  samples  of  a  real  sine  wave  with  initial  phase  4>s 
plus  additive  uncorrelated  noise  samples  e  (n) .  The  angular  frequency  U)s 
is  given  by  8 


u  =  2irf  At  (62) 

s 

where  At  is  the  sampling  interval  (sec)  and  f  is  the  signal  frequency  (Hz) . 
Equation  (61)  can  easily  be  rewritten  in  the  form  of  (52)  by  letting  AQ  =  0.5, 
r=l,  4>0  =  0,  A(J>  =  4>s  -  tt/2,  0)o  =  0  and  Aoi  “  ws.  Then  the  phase  difference 
between  the  two  complex  components  of  the  sine  wave  reckoned  at  the  middle  of 
the  data  record  is,  according  to  (55),  given  by 

A*mid  =  (N_1)a,s  +  2V"  (63) 

The  results  of  Fougere  et  al.  have  been  extended  by  allowing  the 
value  of  r2,  the  ratio  of  the  positive  frequency  to  negative  frequency  signal 
amplitude,  to  range  from  1  to  ■»  in  a  series  of  six  steps.  These  steps  are 
denoted  by  the  letters  A-F,  and  the  relevant  signal  parameters  are  summarized 
in  Table  1.  For  all  steps  the  total  signal  power  Ao(r2+r“2)  was  maintained 
constant  and  equal  to  0.5,  the  power  of  a  real  unit-amplitude  sine  wave.  Also, 
the  values  of  4>o  =  0  and  0)o  =  0  were  maintained  for  all  the  trials.  These 
restrictions  do  not  limit  the  generality  of  the  results  obtained.  It  is  clear 
that  an  arbitrary  phase  rotation  of  the  entire  data-set  will  have  no  effect 
on  its  power  spectrum.  It  is  also  clear  that  since  terms  of  the  form  exp(jmoi0) 
can  be  factored  out  of  the  a(m,M)s  and  B(m,M)s,  the  effect  of  non-zero  w0  is 
simply  to  shift  the  estimated  spectrum  along  the  frequency  axis  (in  a  circular 
or  end-around  fashion)  by  the  amount  <a0.  Finally,  it  should  be  noted  that 
for  each  set  of  cases  examined,  the  same  set  of  noise-data  samples  was  used 
with  all  sets  of  signal  data. 
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TABLE  1 

Test-case  Signal  Parameters 


CASE 

r 

pit) 

r*+r-2 

A 

1 

0 

2 

B 

V2 

0.6 

2.5 

C 

y/lO 

0.98 

10.1 

D 

\/l00 

0.9998 

100.01 

E 

\/l000 

0.999998 

1000.001 

F 

OO 

1 

OO 

5.1  CASE  1 

The  signal  data  for  Case  1  consisted  of  21  samples  of  two  complex  sinu¬ 
soids  with  angular  frequencies  uj  2  °  ±2n/20,  so  that  Aoj  =>  to^.  A<t>micj  was 
stepped  from  -it  to  +tt  in  increments  of  2tt/9  radians,  so  that  in  all  instances 
cosA(J>micj  f  0.  All  spectra  were  estimated  using  (14)  and  a  length  20  (M=19) 
Burg  PEF. 

For  Case  1A  (r=l)  the  situation  corresponds  closely  to  that  of  Fougere 
et  al.'s  Case  2,  where  21  samples  of  a  1  Hz  real  sine  wave  sampled  at  inter¬ 
vals  of  At  =  0.05  s  were  analyzed,  and  <f>8  was  stepped  from  0  to  it  radians  in 
steps  of  tt/9  radian  (20°).  However,  they  used  uncorrelated  real  noise  samples 
uniformly  distributed  over  the  range  [-10“  ,10“^] f  thus  having  a  slightly 
higher  SNR  (82  dB)  than  that  used  here  (77  dB) .  Also,  they  analyzed  their 
data  using  a  real-arithmetic  implementation  of  the  Burg  algorithm. 

The  effects  of  the  relaxation  of  the  conjugate  symmetry  forced  on  the 
spectrum  by  the  use  of  real  data  and  real  arithmetic  analysis  are  apparent 
as  asymmetries  in  Figure  3,  which  shows  a  perspective  projection  of  the  set 
of  10  spectra  viewed  as  though  looking  between  two  parallel  rows  of  pillars. 
A^mid  is  increasing  "into"  the  page,  and  u  is  increasing  from  left  to  right 
as  indicated.  The  vertical  scale  is,  of  course,  distorted  by  the  perspective 
projection;  the  long  and  short  vertical  bars  indicate  20  dB  variation  in 
power  spectral  density  in  the  first  and  tenth  spectra  respectively. 

Figures  4  to  8  show  similarly  displayed  spectra  for  Cases  1B-F,  where 
r  is  increasing  as  shown  in  Table  1.  in  all  cases  the  signal  frequencies 
were  accurately  estimated  by  unsplit  spikes,  although  there  was  some  varia¬ 
bility  of  the  magnitude  of  the  spikes,  particularly  for  the  larger  values  of 
r.  However,  examination  of  the  locations  of  the  zeroes  of  the  PEFs  showed 
that  in  all  cases  the  zeroes  of  the  PEF  varied  by  no  more  than  2x10“^  radians 
from  the  known  locations  of  the  signal  poles. 

The  suggestion  of  Johnson  and  Andersen15  of  computing  and  examining 
the  residues  of  the  poles  of  the  spectral  estimator  (14)  showed  that  the 


Figure  3.  Case  1  A.  Estimated  spectral  power  vs.  angular  frequency  and  A0mj(J.  =  -ojj  =  2rr/20.  N  =  21 
M  =  19.  A0mjd  stepped  from  -n  to  n  in  increments  of  2n/9  radians,  r  =  1.  (Perspective  projection.) 


Figure  4.  Case  IB.  Estimated  spectral  power  vs.  angular  frequency  and  A0mj..  u,  -Wj  2rr/20.  N  21. 
M  19.  N>nwi  steptxHl  from  -it  to  n  in  increments  of  2tr/9  radians,  r  y/2.  (Perspective  projection.) 


Figure  5.  Case  1C.  Estimated  spectral  power  vs.  angular  frequency  and  A#  .j.  cjy  =  -oj*  =  2n/20.  N  =  21. 
M  =  19‘  ^mid  sapped  from  -it  to  n  in  increments  of  2rr/9  radians,  r  *  \J\0.  (Perspective  projection.) 


1  l,JU"‘  6‘  Case  ,D-  Estimated  spectral  power  vs.  angular  frequency  and  w,  =  -w,  =  2rr/20.  N  =  21. 

M  -  19.  A0mjd  stepped  from  -it  to  it  in  increments  of  2w/9  radians,  r  =  VvOO.  (Perspective  projection.) 


Figure  7.  Case  IE.  Estimated  spectral  power  vs.  angular  frequency  and  A0mj(j-  <*>•)  =  -<*>2  =  2rr/20.  N  =  21 
M  =  19.  stepped  from  -n  to  n  in  increments  of  2n/9  radians,  r  =  \/T6bO.  (Perspective  projection.) 


Figure  8.  Case  IF.  Estimated  spectral  power  vs.  angular  frequency  and  A0mjd.  w,  = -u>2  =  2tr/20.  N  -  21 
M  19-  ^mid  st,,PPed  from  -*  to  tr  in  increments  of  2?r/9  radians.  r  =  °°.  (Perspective  projection.) 
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powers  of  the  signal  poles  were  well  represented  by  the  real  parts  of  the 
residues  calculated  at  the  estimated  locations  of  the  signal  poles  (even  when 
the  spectra  showed  gross  variability  in  level,  as  in  Case  ID  (Figure  5)). 
Further,  the  residues  corresponding  to  such  apparently  spurious  spikes  as  are 
located  slightly  to  the  right  of  the  signal  poles  in  Figures  3  to  8  were 
negligible  (always  less  than  10“^  of  the  signal  power). 

The  absence  of  line-splitting  is  fully  consistent  with  the  theory  of 
Section  4.2  as  discussed  in  Section  4.3,  since,  for  these  data,  (58)  is  satis¬ 
fied  in  all  instances:  i.e.,  G(N-l,Au>)  =  0  because  sin[20x(27i /20)  ]  =  0. 
Therefore  at  the  stage  M=2  of  the  Burg  algorithm  8(2,2)  was  estimated  accord¬ 
ing  to  (60)  as  -1,  and  the  error  output  signals  being  processed  at  stages  of 
the  algorithm  beyond  M=2  were  essentially  white  noise  of  low  power.  It 
appears  that  proceeding  beyond  the  stage  M=2  did  not  corrupt  the  estimated 
locations  of  the  signal  poles,  at  least  when  sufficient  dynamic  range  (13 
significant  figures;  floating-point  arithmetic)  is  maintained  in  the  processor. 

Further  examination  of  the  signal  powers  estimated  by  the  real  part  of 
the  residues  of  the  estimated  signal  poles  revealed  that,  while  the  ratio  r2 
of  the  signal  powers  was  accurately  maintained  for  all  cases  examined,  the 
total  signal  power  showed  the  modulation  by  the  factor  cosA4>mi(j,  and  that  the 
depth  of  modulation  was  inversely  proportional  to  (r2+r~^),  as  predicted  by 
(53).  These  observations  provide  further  confirmation  of  the  validity  of  the 
theoretical  results  of  Section  4.2. 


5.2  CASE  2 

The  signal  data  for  Case  2  consisted  of  6  samples  of  two  complex  sinu¬ 
soids  with  angular  frequencies  2  =  ±2tt/4,  so  that  Auj^uj.  A4>mid  was  varied 
from  — 5tt/2  to  — it/ 2  inclusive  in  §1  steps  of  2ir/90  radians  (4°).  All  spectra 
were  estimated  using  (14)  and  a  length  6  (M»5)  PEF.  Thus  Case  2A  (r=l) 
corresponds  to  Case  3  of  Reference  8,  where  unit-amplitude  real  sine  waves  of 
frequency  f  =  5  Hz  were  each  sampled  six  times  at  intervals  At  =  0.05  s  as  the 
initial  phase  4>s  was  stepped  between  0  and  180°  in  increments  of  2°. 

For  this  case  (N-l)Aw  =*  5x  2tt/4  =  5ir/2  so  that  the  data  record  contained 
an  odd  number  of  quarter  cycles  and  G(N-l,Aw)  =  0.2  +  0.  Therefore,  in  light 
of  (56)  and  (60)  it  would  be  expected  that  the  Burg  PEF  would  be  biased  and 
display  dependence  on  cosA^^. 

Substitution  of  N*6  and  Aw*>ir/2  into  (53)  shows  that  1^(0)  is  an  unbiased 
estimate  of  the  input  signal  plus  noise  power,  so  that  fluctuations  in  the 
estimated  signal  powers  made  by  the  residue  method  are  not  expected.  For  this 
case  (56)  reduces  to 


(64) 


(65) 
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Consideration  of  (64)  and  (65)  shows  that  here  the  Burg  PEF  should  have 
greatest  bias  for  r=l,  when  (r2+r~^)  has  its  minimum  value  of  2,  and  that  the 
bias  should  vanish  for  r=",  when  p(r)=l  and  | < 62 ( 1 » 1 ) > |  °  1-4x10“®.  For  the 
case  r  =  ",  (60)  and  hence  (65)  are  not  valid. 


Examination  of  Figures  9  to  14  supports  all  these  conjectures.  These 
figures  show  orthographic  projections  of  the  91  spectra,  with  u>  ranging  from 
-it  to  +7T  across  the  page,  and  A4>mid  increasing  "into"  the  page.  The  labelled 
vertical  bar  to  the  left  of  the  spectra  indicates  a  variation  of  20  dB  in 
power  spectral  density.  It  is  clear  that  there  Is  no  line-splitting  for 
A'hnlil  =  “5ti/2,  —  3  ti  /  2  and  -tt/2,  and  that  the  splitting  shows  a  quas  I -cos  I  hu¬ 
so  Ida  I  dependence  on  as  suggested  by  the  form  of  (64).  The  splitting 

became  less  severe  as  again  as  might  be  inferred  from  (64),  until  for 

Case  2K  the  weaker  signal  pole  was  correctly  estimated  as  a  single  pole. 
Examination  of  the  zeroes  of  the  PEFs  and  the  residue  powers  showed  measurable 
splitting  of  the  stronger  signal  poles  even  for  this  case.  Case  2F  of  course 
showed  no  line  splitting,  since  there  was  then  only  one  signal  pole. 


The  "banding"  effect  visible  most  clearly  in  Figure  9  (Case  2A)  and  to 
a  decreasing  extent  in  subsequent  figures  can  be  explained  on  the  basis  that 
when  cosA4>mid  =  0,  (62(2*2)]  s  1  so  that  the  output  error  power  was  greatly 
reduced  in  those  cases.  This  caused  a  shift  in  the  level  of  the  spectrum,  as 
can  be  seen  from  the  dependence  through  (8b)  of  the  numerator  of  (14)  on  this 
quantity.  This  also  explains  the  obvious  drop  in  spectral  level  in  Figure  14, 
where  1 62 (1*1)  I  ~  1  for  all  values  of  A<J>mid. 


Figurn  11.  Case  2C.  Estimated  spectral  power  vs.  angular  frequency  and  A0mj(j.  w1  -co2  ~  tt/2.  N  6.  M 
^mid  s,t'PfX!d  ,rom  -5»r/2  to  -it/ 2  in  increments  of  2rr/90  radians,  r  \/l0.  (Orthographic  projection.) 


Figure  13.  Case  2E.  Estimated  spectral  power  vs.  angular  frequency  and  A0micj.  wj  =  -u>2  =  tr/2.  N  =  6.  M  =5. 
^mid  *tePPed  from  ~5w/2  to  -n/2  in  increments  of  2tt/90  radians,  r  =  \/l000.  (Orthographic  projection). 
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Figure  14.  Case  2F.  Estimated  spectral  power  vs.  angular  frequency  and  A^j^.  W|  =  -u>2  =  »r/2.  N  =  6.  M  =  5. 
A0m|d  stepped  from  -57t/2  to  -n/2  in  increments  of  2tt/90  radians,  r  =  °°.  (Orthographic  projection.) 


Examination  of  the  residue  powers  showed  that,  when  line-splitting 
occurred,  a  significant  portion  of  the  signal  power  was  accounted  for  by  each 
pole  of  a  split  pair;  in  fact  for  Case  2A  and  ] cosA0mi(i|  «•  1,  70%  of  the 
signal  powers  appeared  at  the  more  severely  deviated  poles,  and  only  30%  at 
less  severely  deviated  poles.  That  the  greater  portion  of  the  signal  power 
was  associated  with  the  more  deviated  pole  appeared  for  these  data  to  be  true 
in  general.  The  present  theory  makes  no  prediction  as  to  why  this  should  be 
the  case,  although  in  principle  the  theory  for  the  infinite  SNR  model  could 
be  extended  to  do  so. 


5.3  CASE  3 

The  signal  data  for  Case  3  consisted  of  25  sets  of  101  samples  of  two 
complex  sinusoids  with  ”  2ir  cases.  Again  wQ“0  was  chosen,  and 

Aw  was  stepped  from  2ir  x  0.0125  to  2tt  x  0.4925  inclusive  in  increments  of 
2v  x  0.02.  Thus  Case  3A  parallels  Case  4  of  Reference  8,  where  101  samples 
were  taken  at  Intervals  At  -  0.01  s  of  real  unit-amplitude  sine  waves  with 
$g  -  */4  and  fa  stepped  from  1.25  Hz  to  49.25  Hz  inclusive  in  steps  of  2  Hz. 
In  all  cases  the  spectra  were  estimated  using  (14)  and  a  length  25  (M-24) 
Burg  PEF. 

Figures  15  to  20  (Cases  3A  to  3F)  show  the  spectral  estimates  obtained 
from  the  data.  These  figures  are  again  orthographic  projections  with  Aw 
increasing  "into"  the  page.  Any  comments  on  the  detailed  structure  of  the 
line-splitting  shown  would  necessarily  be  speculative,  but  some  general 
observations  can  be  made. 


Figure  17.  Case  3C.  Estimate  spectral  power  vs.  angular  frequency.  A$mid  =  2 it.  N  =  101.  M  =  24.  u>i  =  -o 
stepped  from  2 n  x  0.0125  to  2n  x  0.4925  in  increments  of  2ir  x  0.02.  r  =  >/i0.  (Orthographic  projection.) 


Figurn  18.  Case  30.  Estimato  spectral  power  vs.  angular  frequency.  A0mjd  =  2>r.  N=101.  M  =  24.  = -u 

stepped  horn  2tr  x  0.0125  to  2n  x  0.4925  in  increment*  of  2n  x  0.02.  r  =  VlOO.  (Orthographic  projection.) 
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Fiifurc  19.  Case  3E.  Estimate  spectral  power  vs.  angular  frequency.  A0mjd  ‘  2rr.  N  101.  M  24.  u>j  -u>2 
stepped  from  2rr  x  0.0125  to  2tr  x  0.4925  in  increments  of  2»r  x  0.02.  r  =  VlOOO.  (Orthographic  projection.) 


Figure  20.  Case  3F.  Estimate  sp<<ctral  power  vs.  angular  frequency.  2n.  N  101.  M  24.  toj 

stf’PIM'ri  liom  7n  x  0.012b  lo  ?n  x  0.492b  in  iitcimncnls  ol  2 n  x  0.02.  i  (Oi  th<M|i<iphic  projection.) 

The  first  remark  is  that  line-splitting  appears  to  become  less  severe 
as  r  was  increased  in  value,  as  examinations  of  (56)  and  (60)  might  suggest, 
and  in  fact  line-splitting  vanished  for  r*><»  as  discussed  in  Section  5.2. 

The  second  remark  concerns  the  possible  dependence  of  the  spectral 
level  on  the  values  of  G(N-l,Aw)  and  G(N-2,Aw).  These  values  are  given  for 
the  plotted  spectra  in  Table  2.  It  is  interesting  to  note  that  the  minimum 
spectral  level  occurred  at  the  minimum  values  for  |G(N-l,Au) | ,  |G(N-2,Aw)| 
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and  | cosAw |,  and  that  higher  spectral  levels  were  observed  when  G(N-2, Aw) 

<  0,  or  Aw  >  2ir  x  0.25.  Examination  of  (60)  shows  that  as  Aw  exceeds  the 
value  ir/2  certain  terms  change  sign  in  such  a  manner  as  to  decrease  the 
magnitude  of  &2(2,2)  and  thus  perhaps  to  increase  the  magnitude  of  the  numera¬ 
tor  of  (14)  through  the  relation  (8b). 

Finally,  for  Case  3F  it  is  again  observed  that  the  spectral  level  drops 
as  r-*»  and  the  bias  term  in  (56)  vanishes,  similar  to  Case  2F.  For  Case  3F 
only,  the  estimated  poles  accurately  reflected  the  true  signal  pole  locations, 
were  unsplit  and  had  correct  residue  powers.  For  all  other  cases,  examination 
of  the  locations  of  the  poles  of  the  Burg  spectral  estimate  showed  the 
existence  of  multiple  poles  with  significant  residue  power  in  the  vicinity  of 
the  true  locations  of  the  two  signal  poles. 


TABLE  2 

Values  of  (tsno/2tt)  x  100,  G(N-  f.AojJ  and  G(N-2.ku>)  for  Data  from  Case  3 


2m)  x  100 

G(N-1.Aw> 

G(N-2^ir) 

co*Au> 

1.25 

0.1275 

0.1823 

0.9969 

3.25 

0.0493 

0.0488 

0.9792 

5.25 

0.0309 

0.0295 

0.9461 

7.25 

0.0227 

0.0206 

0.8980 

9.25 

0.0182 

0.0125 

0.8358 

11.25 

0.0154 

0.0118 

0.7604 

13.25 

0.0135 

0.0092 

0.6730 

15.25 

0.0122 

0.0071 

0.5750 

17.25 

0.0113 

0.0053 

0.4679 

19.25 

0.0107 

0.0038 

0.3535 

21.25 

0.0103 

0.0024 

0.2334 

23.25 

0.0101 

0.0011 

0.1097 

25.25 

0.0100 

-0.0002 

-0.0157 

27.25 

0.0101 

-0.0014 

-0.1409 

29.25 

0.0104 

-0.0028 

-0.2639 

31.25 

0.0108 

-0.0042 

-0.3827 

33.25 

0.0115 

-00058 

-0.4955 

35.25 

0.0125 

-0.0076 

-0.6004 

37.25 

0.0135 

-0.0098 

-0.6959 

39.25 

0.0160 

-0.0126 

-0.7804 

41.25 

0.0191 

-0.0165 

-0.8526 

43.25 

0.0243 

-0.0224 

-0.9114 

45.25 

0.0340 

-0.0328 

-0.9558 

47.25 

0.0582 

-0.0579 

-0.9551 

49.25 

0.2123 

-0.2142 

-0.9889 
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6.  SUMMARY  AND  CONCLUSIONS 

The  theoretical  properties  of  autoregressive  spectral  analysis  schemes 
have  been  analyzed  when  the  signal  under  investigation  is  known  to  be  compri¬ 
sed  of  one  or  two  complex  sinusoids  in  additive  white  noise.  This  latter  case 
includes  as  a  special  case  data  comprising  a  single  real  sine  wave  in  additive 
white  noise.  It  has  been  shown  that  when  the  autocorrelation  of  the  signal  is 
known,  the  frequency  of  a  single  complex  sinusoid  can  always  be  extracted, 
independent  of  the  signal-to-noise  ratio  (SNR),  provided  enough  samples  of  the 
autocorrelation  are  available.  It  was  also  shown  that  the  Burg  algorithm 
correctly  extracts  the  frequency  of  a  single  sinusoid  in  additive  white  noise 
from  the  complex  amplitude  time-series  data,  provided  the  SNR  is  sufficiently 
high. 


The  situation  when  two  complex  sinusoids  are  present  is  much  more 
complicated.  It  was  found  to  be  fairly  difficult  to  derive  general  equations 
describing  the  spectrum  for  the  known  autocorrelation  (KA)  case,  even  for  a 
two-pole  autoregressive  model.  Nevertheless  these  equations  served  as  a  use¬ 
ful  touchstone  for  the  extremely  complicated  Burg  equations  for  the  analysis 
of  time-series  amplitude  data.  Detailed  theoretical  analysis  showed  that, 
unlike  the  KA  case,  the  Burg  spectral  estimate  is  expected  to  be  sensitive 
both  to  the  number  of  cycles  of  the  difference  frequency  between  the  two 
components  contained  in  the  finite-length  data  record,  and  in  particular  to 
the  relative  phase  difference  between  the  two  complex  sinusoidal  components 
at  the  middle  of  the  data  record.  Finally,  simulation  results  were  shown  to 
be  fully  compatible  with  the  conjectured  basis  of  line-splitting  presented 
here. 
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APPENDIX  A 

A  Bias  Term 

The  bias  term  B2(1,N)  found  in  the  expression  for  <82(1,1)>  is: 
B2(1,N)  -  (N— +  [ (N-2)/ (N-l) ]  X 

[  cos  Aw  +  jp(r)sinAu  +  cosA$midcosAuG(N-2,Au))/ (r2*^  )] 

/[cosAu  +  jp(r)sinAw  +  2cosA4>m^dcosAu>G(N-l,Aw)/(r2+r_S] 


(Al) 


